
source(script.dir%&%'calibrate_model.R') # generated Net.Demand.Adjustment object

# Run simulation with US and demand adjustment
# Also run with 1% oil and gas price shocks to compute effective elasticities
sim.base = excess.demand(p.adj.oil = 1, p.adj.gas=1, 
                         weo.proj=weo.projections, price.path.dummy=price.path,
                         base.spuds.adj = base.spuds.adj,
                         spudregs.use = lm.regs,
                         oil.dem.elast = oil.dem.elast.use,
                         gas.dem.elast = gas.dem.elast.use,
                         Net.Demand.Adjustment = Net.Demand.Adjustment, Demand.Share=0.5, 
                         US.adjustment.factors = US.adj,
                         base.prices = weo.projections[,.(date, wti=wti.base, brent=brent.base, hh=hh.base)],
                         sim.start=sim.start.date, sim.end=sim.end.date,
                         global.gas=global.gas.market)

sim.base.oil.shock = excess.demand(p.adj.oil = 1.01, p.adj.gas=1.01,
                                   weo.proj=weo.projections, price.path.dummy=price.path,
                                   base.spuds.adj = base.spuds.adj,
                                   spudregs.use = lm.regs,
                                   oil.dem.elast = oil.dem.elast.use,
                                   gas.dem.elast = gas.dem.elast.use,
                                   Net.Demand.Adjustment = Net.Demand.Adjustment, Demand.Share=0.5, 
                                   US.adjustment.factors = US.adj,
                                   base.prices = weo.projections[,.(date, wti=wti.base, brent=brent.base, hh=hh.base)],
                                   sim.start=sim.start.date, sim.end=sim.end.date,
                                   global.gas=global.gas.market)

sim.base.gas.shock = excess.demand(p.adj.oil = 1, p.adj.gas=1.01, 
                                   weo.proj=weo.projections, price.path.dummy=price.path,
                                   base.spuds.adj = base.spuds.adj,
                                   spudregs.use = lm.regs,
                                   oil.dem.elast = oil.dem.elast.use,
                                   gas.dem.elast = gas.dem.elast.use,
                                   Net.Demand.Adjustment = Net.Demand.Adjustment, Demand.Share=0.5, 
                                   US.adjustment.factors = US.adj,
                                   base.prices = weo.projections[,.(date, wti=wti.base, brent=brent.base, hh=hh.base)],
                                   sim.start=sim.start.date, sim.end=sim.end.date,
                                   global.gas=global.gas.market)

# Output US oil and gas elasticities for sensitivity case
us.elasts = (log(sim.base.oil.shock$us.sim$OilProduction_AllWells[,-1]/
                   sim.base$us.sim$OilProduction_AllWells[,-1])/log(1.01))[,liq_Total]
us.elasts = data.table(date=sim.base.oil.shock$us.sim$OilProduction_AllWells$date, us.oil.elast=us.elasts)
us.elasts[, us.gas.elast := (log(sim.base.gas.shock$us.sim$GasProduction_AllWells[,-1]/
                                   sim.base$us.sim$GasProduction_AllWells[,-1])/log(1.01))[,gas_Total]]

# fwrite(us.elasts, file=output.dir%&%'/US_avg_elasticities_'%&%prices.to.use%&%'_'%&%dem.elasts.to.use%&%'_'%&%ifelse(us.supp.elast,'ROW-elasts-US','ROW-elasts-IEA')%&%'_'%&%Sys.Date()%&%'.csv')

library(rootSolve)
simulate.policy = function(roy.fed.new.on.sim = NULL,
                           roy.fed.new.off.sim = NULL,
                           onshore.only=FALSE, # only applies to carbon adders and fed leasing. for RRs, specify directly
                           carb.add.fed.sim = NULL, # the single value will override oil/gas specific ones
                           carb.add.fed.sim.oil = 0, 
                           carb.add.fed.sim.gas = 0, 
                           carb.add.nonfed.sim = NULL, # the single value will override oil/gas specific ones
                           carb.add.nonfed.sim.oil = 0, 
                           carb.add.nonfed.sim.gas = 0, 
                           fed.leasing.ban.sim = FALSE, ban.lag.sim=0, # lag in ban taking effect, in years
                           oil.dem.elast = oil.dem.elast.use,
                           gas.dem.elast = gas.dem.elast.use,
                           global.gas.sim = global.gas.market,
                           Net.Demand.Adjustment.sim = Net.Demand.Adjustment,
                           tol = 0.01) {
  if (!is.null(carb.add.fed.sim)) carb.add.fed.sim.oil <- carb.add.fed.sim.gas <- carb.add.fed.sim
  if (!is.null(carb.add.nonfed.sim)) carb.add.nonfed.sim.oil <- carb.add.nonfed.sim.gas <- carb.add.nonfed.sim
  # Important: make sure that the options in the following are IDENTICAL
  st = Sys.time()
  library(rootSolve)
  mr.out = multiroot(f = function(p) {
    excess.demand(p.adj.oil = p[1], p.adj.gas = p[2], 
                  weo.proj=weo.projections, price.path.dummy=price.path,
                  base.spuds.adj = base.spuds.adj,
                  spudregs.use = lm.regs,
                  oil.dem.elast = oil.dem.elast,
                  gas.dem.elast = gas.dem.elast,
                  Net.Demand.Adjustment = Net.Demand.Adjustment.sim, Demand.Share=0.5, 
                  US.adjustment.factors = US.adj,
                  base.prices = weo.projections[,.(date, wti=wti.base, brent=brent.base, hh=hh.base)],
                  sim.start = sim.start.date,
                  sim.end = sim.end.date,
                  roy.fed.new.on.sim = roy.fed.new.on.sim, 
                  roy.fed.new.off.sim = roy.fed.new.off.sim, 
                  onshore.only=onshore.only,
                  carb.add.fed.sim.oil = carb.add.fed.sim.oil,
                  carb.add.fed.sim.gas = carb.add.fed.sim.gas,
                  carb.add.nonfed.sim.oil = carb.add.nonfed.sim.oil,
                  carb.add.nonfed.sim.gas = carb.add.nonfed.sim.gas,
                  fed.leasing.ban = fed.leasing.ban.sim, ban.lag=ban.lag.sim,
                  global.gas=global.gas.sim)$Excess.Demand}, 
    start = c(1, 1))
  print(mr.out)
  ed = Sys.time()
  ed.out = excess.demand(p.adj.oil = mr.out$root[1], p.adj.gas = mr.out$root[2], 
                         weo.proj=weo.projections, price.path.dummy=price.path,
                         base.spuds.adj = base.spuds.adj,
                         spudregs.use = lm.regs,
                         oil.dem.elast = oil.dem.elast,
                         gas.dem.elast = gas.dem.elast,
                         Net.Demand.Adjustment = Net.Demand.Adjustment.sim, Demand.Share=0.5, 
                         US.adjustment.factors = US.adj,
                         base.prices = weo.projections[,.(date, wti=wti.base, brent=brent.base, hh=hh.base)],
                         sim.start = sim.start.date,
                         sim.end = sim.end.date,
                         roy.fed.new.on.sim = roy.fed.new.on.sim, 
                         roy.fed.new.off.sim = roy.fed.new.off.sim, 
                         onshore.only=onshore.only,
                         carb.add.fed.sim.oil = carb.add.fed.sim.oil,
                         carb.add.fed.sim.gas = carb.add.fed.sim.gas,
                         carb.add.nonfed.sim.oil = carb.add.nonfed.sim.oil,
                         carb.add.nonfed.sim.gas = carb.add.nonfed.sim.gas,
                         fed.leasing.ban = fed.leasing.ban.sim, ban.lag=ban.lag.sim,
                         global.gas=global.gas.sim)
  
  if (any(ed.out$ExcessDemand > tol*c(100, 600))) stop('Model failed to converge')
  if (ed.out$global.sim[, any(Oil.Inventories<0)]) warning('Warning: Negative oil inventories')
  if (ed.out$global.sim[, any(Gas.Inventories<0)]) warning('Warning: Negative gas inventories')
  print(difftime(ed, st))
  return(list(sim=ed.out, mr=mr.out))
  
}

# Run policies
st = Sys.time()
sim.base = simulate.policy()
sim.18.75.royalty.rate.on = simulate.policy(roy.fed.new.on.sim=0.1875)
sim.25.royalty.rate.on = simulate.policy(roy.fed.new.on.sim=0.25)
sim.25.royalty.rate = simulate.policy(roy.fed.new.on.sim=0.25, roy.fed.new.off.sim=0.25)
sim.carbon.adder.fed.20 = simulate.policy(carb.add.fed.sim=20)
sim.carbon.adder.fed = simulate.policy(carb.add.fed.sim=50)
sim.ban.lag.0 = simulate.policy(fed.leasing.ban.sim = TRUE, ban.lag=0)
ed = Sys.time()
sim.time = difftime(ed, st)
sim.time

sim.list = list(sim.base=sim.base$sim,
                sim.18.75.royalty.rate.on = sim.18.75.royalty.rate.on$sim,
                sim.25.royalty.rate.on = sim.25.royalty.rate.on$sim,
                sim.25.royalty.rate = sim.25.royalty.rate$sim,
                sim.carbon.adder.fed.20=sim.carbon.adder.fed.20$sim,
                sim.carbon.adder.fed=sim.carbon.adder.fed$sim,
                sim.ban.lag.0=sim.ban.lag.0$sim)

sim.name.list = names(sim.list)

mr.list = list(sim.base=sim.base$mr,
               sim.18.75.royalty.rate.on = sim.18.75.royalty.rate.on$mr,
               sim.25.royalty.rate.on = sim.25.royalty.rate.on$mr,
               sim.25.royalty.rate = sim.25.royalty.rate$mr,
               sim.carbon.adder.fed.20=sim.carbon.adder.fed.20$mr,
               sim.carbon.adder.fed=sim.carbon.adder.fed$mr,
               sim.ban.lag.0=sim.ban.lag.0$mr)
length(sim.list); length(mr.list)


subfolder = ifelse(prices.to.use=='futures','/Futures Prices/', '/High IEA Prices/')
if (dem.elasts.to.use=='Base') subfolder = '/Base elasticity cases'%&%subfolder
if (dem.elasts.to.use=='High') subfolder = '/High elasticity cases'%&%subfolder
if (us.supp.elast) subfolder = '/ROW-US elasticity cases'%&%subfolder
if (low.supp.elast) subfolder = '/Low ROW elasticity cases'%&%subfolder
if (!global.gas.market) subfolder = '/Domestic gas market'%&%subfolder
dir.create(output.dir%&%subfolder, recursive = T)

deltas = mapply(FUN=function(s, m) list(calc.deltas(sim.base, list(sim=s, mr=m))),
                sim.list, mr.list)
names(deltas) = sim.name.list

deltas.summary = sapply(deltas, function(x) list(x$deltas[,.(year, Global.Emis.Tot, Fed.Emis.Tot)]))

co2.results.table = do.call(rbind, lapply(deltas.summary[-1], function(x) {
  out = round(c(x[,Global.Emis.Tot],x[,Fed.Emis.Tot]),2)
  names(out) = rep(x[,year], times=2)
  return(out)
}))

co2.results.table = cbind(co2.results.table, sapply(deltas[-1], function(x) paste0(as.character(round(x$leakage,6)*100),'%')))

price.changes = t(sapply(mr.list[-1], function(x) paste0(as.character(100*round(x[[1]]-1, 6)), '%')))
colnames(price.changes) = c('Oil','Gas')

co2.results.table = cbind(price.changes, co2.results.table)

write.csv(co2.results.table, file=output.dir%&%subfolder%&%'Results_Table_'%&%gsub('/','',subfolder)%&%'_'%&%today()%&%'.csv', row.names = T)


# Sensitivity: how if we adjust the gas emissions rate down to reflect EIA AEO, 
# Gas emissions rate: (117+28.55)/2204.62 = 0.06602045
# how does that change net emissions?
{
  # Save gas emis rate
  gas.emis.rate.backup.fed <- gas.emis.rate.fed
  gas.emis.rate.backup.nonfed <- gas.emis.rate.nonfed
  # gas.emis.rate.fed <- gas.emis.rate.nonfed <-  0.0442001914948831 # high vs low case, EIA AEO 2021
  gas.emis.rate.fed <- gas.emis.rate.nonfed <-  0.0396728492514688 # high vs low case, EIA AEO 2021, don't count petroleum emissions

  deltas.low.gas.emit = mapply(FUN=function(s, m) list(calc.deltas(sim.base, list(sim=s, mr=m))),
                               sim.list, mr.list)
  names(deltas.low.gas.emit) = sim.name.list

  deltas.summary.low.gas.emit = sapply(deltas.low.gas.emit, function(x) list(x$deltas[,.(year, Global.Emis.Tot, Fed.Emis.Tot)]))

  co2.results.table.low.gas.emit = do.call(rbind, lapply(deltas.summary.low.gas.emit[-1], function(x) {
    out = round(c(x[,Global.Emis.Tot],x[,Fed.Emis.Tot]),2)
    names(out) = rep(x[,year], times=2)
    return(out)
  }))

  co2.results.table.low.gas.emit = cbind(co2.results.table.low.gas.emit, sapply(deltas.low.gas.emit[-1], function(x) paste0(as.character(round(x$leakage,6)*100),'%')))

  price.changes.low.gas.emit = t(sapply(mr.list[-1], function(x) paste0(as.character(100*round(x[[1]]-1, 6)), '%')))
  colnames(price.changes.low.gas.emit) = c('Oil','Gas')

  co2.results.table.low.gas.emit = cbind(price.changes.low.gas.emit, co2.results.table.low.gas.emit)

  write.csv(co2.results.table.low.gas.emit, file=output.dir%&%subfolder%&%'Results_Table_LowGasEmitNoPetroleumEmit_'%&%gsub('/','',subfolder)%&%'_'%&%today()%&%'.csv', row.names = T)
  # Return emissions rates to their original values
  gas.emis.rate.fed <- gas.emis.rate.backup.fed
  gas.emis.rate.nonfed <- gas.emis.rate.backup.nonfed
}

gov.revenue(sim.base, agg=T)$tot.rev.summary[year==2035, ]
gov.revenue(sim.carbon.adder.fed, agg=T)$tot.rev.summary[year==2035, ]
gov.revenue(sim.base)$tot.rev.summary

output.sheet.names = c("Baseline", "18.75% RR, Onshore",    "25% RR, Onshore",               
                       "25% RR, On+Off",     "$20 Carbon Adder",            "$50 Carbon Adder",        
                        "Moratorium")

library(openxlsx)
wb = createWorkbook('Revenues')
for (i in 1:length(output.sheet.names)) {
  name.temp = output.sheet.names[i]
  addWorksheet(wb, sheetName=name.temp)
  roy = gov.revenue(sim.list[i], agg=T)$roy.rev.summary
  names(roy)[2:9] = 'roy_'%&%names(roy)[2:9]
  carb = gov.revenue(sim.list[i], agg=T)$carb.rev.summary
  names(carb)[2:9] = 'carb_'%&%names(carb)[2:9]
  
  writeData(wb, sheet=name.temp, cbind(roy, carb[,-'year']))
}
saveWorkbook(wb, output.dir%&%subfolder%&%'/Revenues_'%&%prices.to.use%&%'-prices_'%&%dem.elasts.to.use%&%'-elast_'%&%today()%&%'.xlsx', overwrite = T)
rm(wb)

# Write table:
sim.name.list
sims.for.table = 1:length(sim.name.list)
sim.name.list[sims.for.table]

sum.table = mapply(FUN=extract.vars.for.table, d=deltas[sims.for.table], nm=sims.for.table)
sum.table = t(sum.table)
sum.table
sum.table.num = sum.table
write.csv(sum.table, file=output.dir%&%subfolder%&%'Summary_table'%&%'_'%&%today()%&%'.csv')

sum.table[,-c(1,2,7,8)] = round(sum.table[,-c(1,2,7,8)], 0)
sum.table[,c(1,2,7)] = scales::percent(sum.table[,c(1,2,7)])
sum.table[,8] = scales::dollar(as.numeric(sum.table[,8]))
sum.table = sum.table[-1,] # drop baseline

rownames(sum.table) = c('18.75% Royalty Rate (Onshore only)',
                        '25% Royalty Rate (Onshore only)',
                        '25% Royalty Rate (Onshore and Offshore)',
                        '$20 Carbon Adder (2%)',
                        '$50 Carbon Adder (2%)',
                        'Moratorium')
library(stargazer)

# stargazer(sum.table, out = output.dir%&%subfolder%&%'Summary_table_low-ROW-elast'%&%today()%&%'.tex')
stargazer(sum.table, out = output.dir%&%subfolder%&%'Summary_table_'%&%today()%&%'.tex')

sim.name.list[sims.for.table]
scen.labels = c('18.75% Royalty Rate (Onshore only)',
                '25% Royalty Rate (Onshore only)',
                '25% Royalty Rate (Onshore and Offshore)',
                '$20 Carbon Adder (2%)',
                '$50 Carbon Adder (2%)',
                'Moratorium')

# Plot price projections
# dev.off()
pdf(output.dir%&%subfolder%&%'/Price_projections_'%&%prices.to.use%&%'_'%&%today()%&%'.pdf', family='serif', width=11, height=8.5)
par(mai=par('mai')*c(1,1,1,2))
plot(wti.base ~ date, data=weo.projections, type='l', col='black', lwd=3,
     xlab='',ylab='Oil Price ($/barrel)', ylim=c(0,ifelse(prices.to.use=='futures',100, 110)), cex.lab=1.5, cex.axis=1.5)
grid(nx=NA, ny=NULL)
abline(v=seq.Date(as.Date('2020-01-01'),as.Date('2050-01-01'), by='5 years'), col='lightgray', lty=3)
if (prices.to.use=='futures') {
  abline(v=wti.fut[,min(date)], lty=1, col='gray20')
  text(x=wti.fut[,min(date)], labels = 'Historical\nPrices', y=95, pos = 2, cex = 1.25)
  text(x=wti.fut[,min(date)], labels = 'Futures\nPrices', y=95, pos = 4, cex = 1.25)
  abline(v=wti.fut[,max(date)], lty=1, col='gray20')
  text(x=wti.fut[,max(date)], labels = 'Projected\nFutures\nPrices', y=92, pos = 4, cex = 1.25)
}
par(new=TRUE)
plot(hh.base ~ date, data=weo.projections, lty=1, col='blue', 
     lwd=3, type='l', yaxt='n', xaxt='n', ylim=c(0,ifelse(prices.to.use=='futures',100, 110))/10,
     ylab='', xlab='')
axis(4)
mtext('Gas Price ($/mmbtu)', side=4, line=2, cex = 1.5)
legend(ifelse(prices.to.use=='futures','topright','top'), legend=c('WTI (left)','Henry Hub (right)'),
       lty=1, lwd=3, col=c('black','blue'), cex=1.5, bg='white',
       y.intersp=0.75)
dev.off()

x = calc.deltas.sep.leakage(sb=sim.base, sp=sim.ban.lag.0)
x$deltas[.N, .(ROW.Prod.Oil/Fed.Oil, ROW.Prod.Oil/US.Oil)]
x$deltas[.N, .((Nonfed.Oil.On+Nonfed.Oil.Off)/Fed.Oil, ROW.Prod.Oil/Fed.Oil, 
               Global.Prod.Oil/Fed.Oil-1, ROW.Prod.Oil/US.Oil)]

x$deltas[.N, .(US.Oil, ROW.Prod.Oil, US.Oil+ROW.Prod.Oil, Global.Prod.Oil, -US.Oil/ROW.Prod.Oil)]
x$deltas[.N, .(US.Gas, ROW.Prod.Gas, US.Gas+ROW.Prod.Gas, Global.Prod.Gas, -US.Gas/ROW.Prod.Gas)]
req.ratio = x$deltas[.N, -US.Gas/ROW.Prod.Gas]
req.ratio
(req.ratio*gas.emis.rate.fed - 117/2204.62)/(28.55/2204.62) # meaning 4-7x US leakage rates
(req.ratio*gas.emis.rate.fed - 117/2204.62)/(28.55/2204.62)*0.0236 # meaning >=8.5-15% leakage

deltas.sep = mapply(FUN=function(s, m) list(calc.deltas.sep.leakage(sim.base, list(sim=s, mr=m), yrs=2020:2050)),
                    sim.list, mr.list)

deltas.sep$sim.ban.lag.0$deltas[.N,  c(-Global.Emis.Tot, ROW.Emis.Gas, -Global.Emis.Tot+ROW.Emis.Gas)] # net reductions exceed ROW gas emis increase by 136-189 MT/y
# E.g., with -85 in net reductions (accounting for the +51) alongside 51 in increases
# we would need ANOTHER 85 in gas emissions, which is an increase by a factor of 1 + 85/51
# This means ROW gas emissions would have to be larger by a factor of 
req.ratio = deltas.sep$sim.ban.lag.0$deltas[.N,  c(1 + -Global.Emis.Tot/ROW.Emis.Gas)] # ROW gas emissions intensity would have to be 2.67-4.51x larger than US to net to zero.
# req.ratio = deltas.sep$sim.ban.lag.0$deltas[.N,  c(1 + -Global.Emis.Gas/ROW.Emis.Gas)] # ROW gas emissions intensity would have to be 2.67-4.51x larger than US to net to zero.
req.ratio*gas.emis.rate.fed # or 0.17-0.29 tCO2e/mcf 
# Since the carbon content is fixed, to get to this level of emissions intensity by just changing the 
# fugitive methane leakage rate would require 10-19 times larger than the US value (2.36%)
(req.ratio*gas.emis.rate.fed - 117/2204.62)/(28.55/2204.62)
(req.ratio*gas.emis.rate.fed - 117/2204.62)/(28.55/2204.62)*0.023 # meaning >=22%-44% leakage
# Even 22% leakage is more than double any study in the literature. Typically global leakage estimates 
# are in the range of 2-3%.

# Same calc for oil
req.ratio = deltas.sep$sim.ban.lag.0$deltas[.N,  c(1 + -Global.Emis.Tot/ROW.Emis.Oil)]
# req.ratio.oil = deltas.sep$sim.18.75.royalty.rate.on$deltas[.N,  c(1 + -Global.Emis.Tot/ROW.Emis.Oil)]
# req.ratio.oil = deltas.sep$sim.ban.lag.0$deltas[.N,  c(1 + -Global.Emis.Oil/ROW.Emis.Oil)]
req.ratio*oil.emis.rate.fed # 1.15-1.94 tCO2e/barrel. Way above anything realistic.
# Even the dirtiest oil (Canadian tarsands) have all-in emissions of 0.736 tCO2e/barrel 
# https://oci.carnegieendowment.org/

deltas.sep$sim.carbon.adder.fed$deltas[, .(year, Global.Emis.Oil, Global.Emis.Gas, sum(Global.Emis.Oil)/sum(Global.Emis.Gas))]
deltas.sep$sim.ban.lag.0$deltas[, .(year, Global.Emis.Oil, Global.Emis.Gas, sum(Global.Emis.Oil)/sum(Global.Emis.Gas))]
# RFF green
rff.b = rgb(4, 39, 60, maxColorValue=255) # RFF black

for (i in 1:5) {
  if (i==1) {d.temp = deltas.sep$sim.18.75.royalty.rate.on$deltas[-c(1,.N)]; scen.nm = '1875RR-On'}
  if (i==2) {d.temp = deltas.sep$sim.25.royalty.rate$deltas[-c(1,.N)]; scen.nm = '25RR'}
  if (i==3) {d.temp = deltas.sep$sim.carbon.adder.fed$deltas[-c(1,.N)]; scen.nm = 'CarbonAdder50'}
  if (i==4) {d.temp = deltas.sep$sim.ban.lag.0$deltas[-c(1,.N)]; scen.nm = 'Moratorium'}
  if (i==5) {d.temp = deltas.sep$sim.carbon.adder.fed.20$deltas[-c(1,.N)]; scen.nm = 'CarbonAdder20'}
  
  # dev.off()
  pdf(output.dir%&%subfolder%&%'/ProdChangeBreakdown-'%&%scen.nm%&%'_'%&%dem.elasts.to.use%&%'-elasts_'%&%prices.to.use%&%'-prices_'%&%today()%&%'.pdf',
      family='serif', width=11, height=10)
  par(mfrow=c(2,1))
  den = 30
  par(mar=c(1,1,1,3)*par('mar'))
  bar = barplot(height=t(as.matrix(d.temp[,.(Fed.Oil.On, Fed.Oil.Off, Nonfed.Oil.On, Nonfed.Oil.Off, ROW.Prod.Oil)])), 
                beside=T, names.arg=d.temp$year, col=c(my.cols[c(2,2,1,1)], rff.b), density=c(NA,den,NA,den,NA), angle=30,
                lwd=2, xlab='Year', ylim=c(-1.5,1.5), cex.main=1.5, cex.axis=1.25, cex.lab=1.5, cex.names = 1.5,
                ylab='Change in Oil Production (mb/d)', main='Change in Oil Production, by Region')
  grid()
  barplot(height=t(as.matrix(d.temp[,.(Fed.Oil.On, Fed.Oil.Off, Nonfed.Oil.On, Nonfed.Oil.Off, ROW.Prod.Oil)])), 
          ylab='', xlab='', xaxt='n', yaxt='n',
          beside=T, names.arg=d.temp$year, col=c(my.cols[c(2,2,1,1)], rff.b), density=c(NA,den,NA,den,NA), angle=30, add=T)
  par('usr')[1:2]
  lines(d.temp$Global.Prod.Oil ~ bar[3,], data=d.temp, lwd=3, col=3)
  
  legend('bottomleft', legend=c('Federal US Onshore','Federal US Offshore',
                                'Nonfederal US Onshore','Nonfederal US Offshore', 'ROW'),
         fill=c(my.cols[c(2,2,1,1)], rff.b), density=c(NA,den,NA,den,NA), angle=30, bg='white')
  legend('topleft', lty=1, lwd=3, col=3, legend='Net Global', bg='white')
  yrg = par('usr')[3:4]
  yrg = seq(yrg[1], yrg[length(yrg)], length.out=7)
  yrg.2 = seq(yrg[1], yrg[length(yrg)], length.out=7)*oil.emis.rate.fed*365.25
  axis(4, at=yrg, labels = round(yrg.2, 0), cex.lab=1.5, cex.axis=1.25)
  axis(2, at=yrg[2], labels = yrg[2], cex.lab=1.5, cex.axis=1.25)
  axis(4, at=yrg[2], labels = round(yrg.2[2], 0), cex.lab=1.5, cex.axis=1.25)
  mtext('Change in Emissions\n(MMTCO2e/year)', side=4, cex=1.5, line=4)
  
  # Plot gas
  bar = barplot(height=t(as.matrix(d.temp[,.(Fed.Gas.On, Fed.Gas.Off, Nonfed.Gas.On, Nonfed.Gas.Off, ROW.Prod.Gas)])), 
                beside=T, names.arg=d.temp$year, col=c(my.cols[c(2,2,1,1)], rff.b), density=c(NA,den,NA,den,NA), angle=30,
                lwd=2, xlab='Year', ylim=c(-8, 6), cex.main=1.5, cex.axis=1.25, cex.lab=1.5, cex.names = 1.5,
                ylab='Change in Gas Production (bcf/d)', main='Change in Gas Production, by Region')
  grid()
  barplot(height=t(as.matrix(d.temp[,.(Fed.Gas.On, Fed.Gas.Off, Nonfed.Gas.On, Nonfed.Gas.Off, ROW.Prod.Gas)])), 
          ylab='', xlab='', xaxt='n', yaxt='n',
          beside=T, names.arg=d.temp$year, col=c(my.cols[c(2,2,1,1)], rff.b), density=c(NA,den,NA,den,NA), angle=30, add=T)
  par('usr')[1:2]
  lines(d.temp$Global.Prod.Gas ~ bar[3,], data=d.temp, lwd=3, col=3)
  
  
  yrg = par('usr')[3:4]
  yrg = seq(yrg[1], yrg[length(yrg)], length.out=8)
  yrg.2 = seq(yrg[1], yrg[length(yrg)], length.out=8)*gas.emis.rate.fed*365.25
  axis(4, at=yrg, labels = round(yrg.2, 0), cex.lab=1.5, cex.axis=1.25)
  axis(4, at=yrg[c(2,4)], labels = round(yrg.2[c(2,4)], 0), cex.lab=1.5, cex.axis=1.25)
  mtext('Change in Emissions\n(MMTCO2e/year)', side=4, cex=1.5, line=4)
  dev.off()
  mean(d.temp$Global.Prod.Oil)*365.25*oil.emis.rate.fed
  mean(d.temp$Global.Prod.Gas)*365.25*gas.emis.rate.fed
}

# Graph % federal reductions by policy and year
fed.red.by.year = mapply(FUN=function(s, m) {
  x = calc.deltas(sim.base, list(sim=s, mr=m), yrs = 2020:2050)$policy[, (Fed.Emis.Tot)]/
    calc.deltas(sim.base, list(sim=s, mr=m), yrs = 2020:2050)$base[, (Fed.Emis.Tot)]-1
  return(x)
}, sim.list, mr.list)
fed.red.by.year = data.table(fed.red.by.year)
fed.red.by.year = cbind(year=2020:2050, fed.red.by.year)

# Check:
fed.red.by.year[, sim.ban.lag.0]/
  (calc.deltas(sim.base, sim.ban.lag.0, yrs = 2020:2050)$policy[, (Fed.Emis.Tot)]/calc.deltas(sim.base, sim.ban.lag.0, yrs = 2020:2050)$base[, (Fed.Emis.Tot)]-1)

sim.name.list[sims.for.table]

my.cols.range = c('#04273c','#74645e','#ff6663','#765ea5','#88c4f4','#50b161')
names(my.cols.range) = sim.name.list[sims.for.table[-c(1,7)]]
pdf(output.dir%&%subfolder%&%'/Federal_emissions_reductions_percent_'%&%today()%&%'.pdf', family='serif', width=14, height=8)
plot(sim.18.75.royalty.rate.on*100 ~ year, data=fed.red.by.year, 
     ylim=c(-100,0), xlim=c(2020, 2055), xaxt='n', type='n', lwd=2, col=my.cols.range[1],
     xlab='Year',ylab='Percent Reduction in Annual Federal Oil and Gas Emissions',
     cex.axis=1.5, cex.lab=1.5)
axis(1, at=seq(2020, 2050, by=5), cex.axis=1.5)
grid()
lines(sim.18.75.royalty.rate.on*100 ~ year, data=fed.red.by.year, ylim=c(-100,0), lwd=3, col=my.cols.range[1])
lines(sim.25.royalty.rate.on*100 ~ year, data=fed.red.by.year, ylim=c(-100,0), lwd=3, col=my.cols.range[2])
lines(sim.25.royalty.rate*100 ~ year, data=fed.red.by.year, ylim=c(-100,0), lwd=3, col=my.cols.range[3])
lines(sim.carbon.adder.fed.20*100 ~ year, data=fed.red.by.year, ylim=c(-100,0), lwd=3, col=my.cols.range[4])
lines(sim.carbon.adder.fed*100 ~ year, data=fed.red.by.year, ylim=c(-100,0), lwd=3, col=my.cols.range[5])
lines(sim.ban.lag.0*100 ~ year, data=fed.red.by.year, ylim=c(-100,0), lwd=3, col=my.cols.range[6])
text(x=2050, y=100*as.numeric(fed.red.by.year[year==2050, 
                                              c('sim.18.75.royalty.rate.on',
                                                'sim.25.royalty.rate.on',
                                                'sim.25.royalty.rate',
                                                'sim.carbon.adder.fed.20',
                                                'sim.carbon.adder.fed',
                                                'sim.ban.lag.0')])+c(0,1,-3.5,0,0,0),
     labels = c('18.75% Onshore RR',
                '25% RR, Onshore Only',
                '25% RR, Onshore &\nOffshore',
                '$20 Carbon Adder (2%)',
                '$50 Carbon Adder (2%)',
                'Leasing Ban'), col=my.cols.range, pos = 4, cex = 1.2)
dev.off()
